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ABSTRACT 


A  finite  difference  model  for  predicting  the  nearshore 
circulation  due  to  wind  and  waves  is  presented  which  attempts  to 
solve  the  same  problem  as  an  earlier  model  created  by  Birkemeier 
and  Dalrymple  (197S) .  Their  model  iteratively  solved  the  linear 
set  of  conservation  equations  of  both  mass  and  momentum,  which 
were  time  averaged  (over  one  wave  period)  and  depth  integrated, 
for  mean  velocities  and  free  surface  displacements.  The  wave 
characteristics  used  in  the  momentum  equations  were  found  using 
the  wave  refraction  and  shoaling  routines,  including  wave-current 
interaction,  developed  by  Noda,  et  al.  (1974).  The  model  also 
included  a  linear  bottom  friction  formulation  as  well  as  a  surface 
wind  stress  capability. 

The  present  model  discussed  herein  includes  the  addition  of 
convective  accelerations,  horizontal  mixing  and  a  quadratic 
bottom  friction  term  in  the  conservation  of  momentum  equations. 
This  bottom  friction  term  is  "exact"  in  the  sense  that  it  includes 
the  velocity  vectors  due  to  both  mean  and  wave-induced  currents. 


The  model  is  applied  to  the  cases  of  a  single  wave  train 
impinging  on  a  plane  beach,  a  barred  profile,  and  a  bottom  with  a 

periodically  spaced  rip  channel.  It  is  also  applied,  in  a 

pr. - _____ - - -• 

simplified  form,  to  the  case  of  two  intersecting  wave  trains  at 
oblique  angles  to  a  plane  beach.  Results  indicate  that  these 
additions  to  the  model  are  important  in  attempts  to  model  the 
circulation  patterns  over  bottom  bathymetries  found  in  nature. 


v 


CHAPTER  I 
INTRODUCTION 


In  recent  years  there  has  evolved  a  special  interest  in  the 
environment  around  our  extensive  coastline.  This  interest  is  a  result 
of  the  growing  usage  of  these  areas  for  industry  and  recreation. 

Those  industries  which  rely  on  water  as  a  primary  mode  of  transporta¬ 
tion  for  both  raw  materials  and  finished  products  are  spending  large 
amounts  of  time  and  money  designing  efficient  access  and  docking 
facilities.  Coastal  communities  are  spending  more  money  to  maintain 
the  condition  of  their  recreational  facilities  located  around 
water  bodies.  As  people  grow  to  depend  on  these  regions  for  their 
livelihood,  they  become  increasingly  concerned  with  changes  to  the 
coastline  and  with  efforts  to  maintain  it  in  its  present  condition. 

The  problem  of  sediment  transport  is  a  primary  concern  of 
industry  and  resort  communities.  Its  effects  are  clearly  visible, 
yet  an  understanding  of  its  causes  and  the  ability  to  predict  it 
accurately  are  just  starting  to  unfold.  An  integral  part  of  the 
sediment  transport  problem  is  an  ability  to  describe  the  flows  which 
serve  to  move  sediment,  such  as  the  creation  of  a  velocity  field  in 
the  surf  zone  due  to  breaking  waves. 
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Noda,  et  al.  (1974)  developed  a  steady  state  model  which 
predicted  the  nearshore  circulation  due  to  waves  and  wind.  Liu  and 
Lennon  (1978)  also  developed  a  steady  state  model  using  finite 
element  techniques.  Birkemeier  and  Dalrymple  (1975),  building  on  the 
efforts  of  Noda,  et  al . ,  created  a  time  dependent  model  to  describe 
nearshore  circulation  due  to  the  same  forces.  The  effort  reflected 
in  this  thesis  is  an  attempt  to  extend  the  work  of  Birkemeier  and 
Dalrymple  one  step  further;  to  develop  a  more  complete  and  accurate 
numerical  model  to  predict  flows  along  a  coastline  by  including  the 
effects  of  convective  accelerations  and  lateral  mixing,  and  by  the 
use  of  a  bottom  friction  term  which  includes  the  velocities  due  to 
both  waves  and  mean  currents. 

There  are  many  advantages  to  using  a  numerical  model  to 
investigate  this  problem.  First  of  all,  the  problem  at  hand  is  so 
complex  that  only  through  numerical  methods  can  it  be  solved. 

Secondly,  the  model  enables  the  user  to  study  an  unsealed  prototype 
situation.  Thirdly,  once  the  model  has  been  established  its  generality 
enables  it  to  be  used  theoretically  on  any  stretch  of  beach.  Finally, 
the  model  arrives  at  solutions  rapidly  when  adapted  to  use  with  a 
high  speed  computer. 

In  order  to  formulate  the  problem  so  that  it  may  be  solved 
numerically  the  following  requirements  must  be  met: 
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(1)  a  set  of  governing  equations  must  be  selected 
to  accurately  describe  the  physical  processes 
at  hand;  the  model  is  only  as  good  as  its 
governing  equations, 

(2)  boundary  conditions  must  be  established  over 
the  region  of  interest,  and 

(3)  initial  conditions  must  be  defined. 


Only  when  these  requirements  are  met  can  we  hope  to  obtain  an  accurate 
solution. 


CHAPTER  II 
GOVERNING  EQUATIONS 


INTRODUCTION 

The  first  step  in  problem  formulation,  as  just  mentioned, 
is  the  selection  of  governing  equations  that  accurately  describe  the 
physical  processes  at  work.  This  numerical  model  uses  as  its 
basis  the  equations  describing  conservation  of  both  mass  and  momentum 
in  a  time  averaged  (over  one  wave  period)  and  depth  integrated 
form. 


The  purpose  of  the  time  averaging  is  to  remove  the  fluctuations 
in  time  due  to  waves.  The  model  deals  only  in  mean  quantities.  The 
reasons  for  depth  integrating  the  equations  are,  again,  to  deal  with 
mean  quantities,  thus  time  over  depth,  and  to  reduce  the  problem 
from  two  horizontal  and  one  vertical  dimension  into  only  two 
horizontal  dimensions.  In  the  process  of  depth  integration  the 
Leibnitz  Rule  was  used  to  remove  partial  derivatives  from  within  an 
integral.  It  is  given  by 


r8  (y) 

a(y) 


3f (x,y) 
3y 


rfil  y) 
'a(y) 


f (x,y)dx  -  f (6 (y) ,y) 


36  (y) 

3y 


f  (a(y)  ,y)  3«  (y) 
3y 


The  remainder  of  this  chapter  is  devoted  to  the  derivation  of 


the  integrated  conservation  of  mass  and  momentum  equations  and  the 
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"forcing"  terms  found  in  the  momentum  equations  themselves.  The 
three  basic  equations  discussed  above  will  be  derived  in  detail. 

The  equations  which  are  solved  to  yield  the  wave  characteristics, 
the  radiation  stresses,  and  the  wind  stress  will  be  discussed 
briefly.  For  a  more  detailed  derivation  the  reader  is  referred  to 
the  work  of  Noda,  et  al.  and  Birkemeier  and  Dalrymple.  The 
bottom  stress  derivation  and  the  formulation  of  the  lateral  mixing 
terms,  however,  will  be  discussed  more  fully  as  they  are  new 
additions  to  the  model. 

BOUNDARY  CONDITIONS 

Certain  kinematic  boundary  conditions  are  used  in  the  formu¬ 
lation  of  the  mass  and  momentum  equations.  The  boundary  conditions 
state  that  if  we  move  with  a  surface  given  by 

F(x,y,z,t)  =  0  , 

then  a  water  particle  cannot  flow  across  the  surface,  otherwise 
the  surface  would  cease  to  exist.  Mathematically,  this  condition 
is  expressed  by  the  total  time  rate  of  change  of  the  function 
F(x,y,z,t)  equals  zero. 


D 

Dt 


(F(x,y,z,t) ) 


=  0 


At  the  free  surface  the  boundary  is  given  by 
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F^x^Zjt)  =  z  -  n(x,y,t)  =  0 

and  at  the  bottom 

F2(x,y,z,tji  =  z  +  h(x,y,t)  =  0 

Therefore,  for  the  kinematic  free  surface  boundary  condition 

(KFSBC) , 


3F  3F  ^  3F  3F 

3t  U  3x  3y  3z 


=  0 


or 


w 


3n  ,  3n  ,  3n 

-  +  u  — L  +  v  - — 

3t  n  3x  n  3y 


(1) 


For  the  bottom  boundary  condition  (BBC)  it  can  be  shown  that, 


in 

3t 


+  u 


l-h  3x  +  v-h  3y 


0 


(2) 


where  u,v,w  are  the  velocity  components  in  the  x,y,  and  z  directions 
and  the  subscripts  denote  the  location  of  a  specific  term  whether  it 
be  at  the  bottom,  z  =  -h,  or  at  the  free  surface,  z  *  n. 


CONTINUITY  EQUATION 

The  general  form  of  the  three  dimensional  continuity  equation 
which  states  that  mass  is  conserved  is. 


j)£_  3(pu)  3(pv)  3(pw) 

3t  3x  3y  3z  u 
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where  p  is  the  density  of  the  fluid. 

Integrating  over  depth  from  z  ■  -h(x,y,t)  to  z  »  n(x,y,t)  and  using 
the  Leibnitz  Rule,  Equation  (3)  becomes 


_± 

3t 


dz  -  p 

n 


3(-h)  _3  fn 

3t  3x 


pudz 


3n  .  3(-h)  . 

-%un  ^  +  p-hu-h  “3T  + 


57  j> 


3n  3  -h) 

-p  v  —  +  p  .v  .  — 7 —  +  p  w  -  p  ,  w  ,  =  ° 

n  n  9y  -h  -h  3y  n  n  ~h  -h 


Simplifying  by  substituting  in  both  the  KFSBC  and  the  BBC, 
the  continuity  equation  can  be  written 


pdz  + 


pudz  + 


pvdz  =  0 


(4) 


Let  both  u  and  v  be  comprised  of  a  time  independent  mean 
flow,  a  wave  induced  flow,  and  an  arbitrarily  fluctuating  component 
(turbulence)  such  that 

u  =  U  +  0  +  u’ 

v=V+v+v' 

It  is  important  to  note  that  the  turbulent  fluctuations  u' 
and  v'  are  of  very  high  frequency  and  by  definition  their  time 
averages  (over  a  short  period  of  tiue)  are  identically  zero. 


1  f 

£  u'dt  =  0 

T  Jo 
1  fT 

£  v'dt  =  0 

T  Jo 


By  substituting  the  above  expressions  for  u  and  v  into 
Equation  (4) ,  time  averaging  over  one  wave  period  so  as  to  eliminate 
the  wave  induced  fluctuations,  and  using  the  definitions  for  the 
time  average  of  the  turbulent  components ,  the  continuity  equation 
can  be  written  as 


g~  {p(h  +  n) }  +  (pU(h  +  h) }  + 

— v 


*  {pV(h  +  D  >  +  37  j 


pvdz  =  0 


where  the  symbol  " - "  denotes  the  time  average  over  one  wave 

period  and  n  the  time  independent  mean  free  surface  displacement. 
Note  that  the  time  average  of  the  vertically  integrated  wave  induced 

A  ^ 

velocities  u  and  v  is  not  zero. 


Defining 


U  s  u  +  u 


V  =  V  +  v 


where  u 


n  - 

pu< 

'  — h 


p(h  +  n) 


r 

pvc 

>  -h 


p  (h  +  n) 
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are  the  mass  transport  velocities  and  substituting  the  total  depth 
D  for  (h  +  n) /  the  time  averaged  depth  integrated  continuity  equation 
is,  in  its  final  form 

I? + -£ *  4 <vB  ■ 0  <6> 

Also  inherent  in  the  derivation  are  the  assumptions  that  the  bottom 
is  constant  with  time  and  the  density  is  constant  in  space  and 
time. 


MOMENTUM  EQUATIONS 

The  x  and  y  momentum  equations  are  manipulated  in  the  same 
way  as  the  continuity  equation  in  order  to  achieve  equations  which 
are  independent  of  wave  induced  oscillations,  i.e.,  they  are 
integrated  over  depth  and  time  averaged  over  a  wave  period.  The 
general  form  of  the  horizontal  momentum  equations  are;  in  the  x 
direction. 


3u  3u2  3uv  3uw  _  1_  3P_  1_  .9Txx  9Tyx  9tzx- 

3t  +  3x  3y  3z  p  3x  p  3x  3y  3z 


(7) 


and  in  the  y  direction, 


3v  3uv  3v2  3vw  1  3P  1  , 9Txy  .  9  Tyy  .  9  Tzy, 

3t  IhT  37~  9z  p"  3y  p  {^T  + 


(8) 


v 

The  x  momentum  equation  wi1 1  be  dealt  with  first.  Integrating 


the  left  hand  side  over  depth  and  using  the  Leibnitz  Rule  to  remove 
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derivatives  from  within  the  integrals  give  term  by  term, 


(ub>  r  it az  *  ai  r uas  ■  %  it  *  u-t> 


P  3u2  .  3  P  2 

t —  dz  =  —  u  dz  -  u 

J-h  3x  J.  n 


2  3n  2  3(-h) 

3x  U-h  3x 


P  3uv  .  3  P  ,  3n 

t —  dz  =  t—  uvdz  -  u  v  t —  +  u  ,v 
J  .  3y  3y  J  .  n  n  3y  -h  -! 


3(-h) 
h  3y 


P  3uw 

J-h  32 


3uw  . 

- —  dz  =  u  w  -  u  ,  w  , 

3z  n  n  ~h  -h 


Rearranging  terms  the  (LHS)  becomes, 


3t 


r udz + » r  “2az * #  £ 


uvdz 


,3n  9n  3n  , 

-u  (77  +  u  -r±  +  v  t - w  ): 

n  3t  n  3x  n  3y  n 


,3h  A  3h  3h  A  , 

_U-h  ITt  +  u-h  37  +  V-h  37  +  W-h} 


which  after  substituting  the  KFSBC  and  the  BBC  simplifies  to. 


(LHS) 


3t 


rn  a  p  5  3  p 

J-h  UdZ  +  *  J_h  U  dZ  +  F  J.h  uvdz  * 


Again  letting  the  total  velocities  u  and  v  be  defined  as 


u  =  U  +  u  +  u' 


v 


V  +  V  +  v' 


Substituting  and  averaging  over  a  wave  period,  the  (LHS)  becomes 
term  by  term. 


if  „  3  f  3  f  .  A  3  f  .. 

*  L  udz  =  *  L  °dz  +  *  L  udz  +  ^  L  u  dz 


-r5-  f  u2dz  «  fV2 dz  +  P  G2dz  +  ■—  P  u,2dz 

3x  J_h  3x  J-h  3x  J. h  3x  _h 


[  2Uudz  +  -r^-  [  2uu'dz  +  (  2{iu'dz 

3x  J-h  3x  J-h  3X  J-h 


J_hUVdZ  =  9?  LUVdZ  +  37  l.ydZ  +  37  J_hUV’dZ 


3  fn  — 3  fn  .»  3  f 

+  t—  vS/dz  +  t—  uvdz  +  t—  uv'dz 

3y  J.h  3y  J.h  3y  J_h 


3y  Lhu'Vdz  +  4  L  u'^dz  +  F  L  u’v'dz 
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Since  the  frequency  of  the  fluctuating  components  is  much  greater 
than  the  wave  frequency,  the  wave  induced  velocity  is  essentially 
constant  relative  to  the  turbulent  fluctuations.  Therefore,  integrals 
which  involve  products  of  a  turbulent  component  and  a  wave  induced 
component  are  zero  in  the  time  average.  Using  this  result  the  (LHS) 
becomes. 


£  11  "az '  * 11"* +  &  11 iaz 


_a 

P  )  3  f 

n  -2 

a 

fn  -2  a  f 

3x 

"  az  ’  17  J 

U  dz  + 

ax 

L  u  dz  +  ^ 

•'-h  J 

-h 

1  -h  ' 

a  rn 

+  ta  J_h20udz 


0  rn  i  0  rn  A  ^  rn 

17  J.h  uvaz  '  17  J.h  mdz  *  17  J_h  Uvdz  *  17  J_h  “vdz 


a  f1  3  f1 

r—  uvdz  +  T—  \1'' 

3y  i-h  3y  Lh 


The  right  hand  side  after  integrating  over  depth  is 


(RHS) 


.if  !3n,„i r 

p  J-h  3x  p  ^-h  3x  p  i-h  3y  p  i-i 


rn 


\ 

! 
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assuming  the  density  is  constant  in  space  and  time.  If  we  also  assume 
an  inviscid  fluid  such  that  no  horizontal  viscous  stress  exists,  then 
and  become  zero.  Integrating  over  a  wave  period  and  invoking 
the  Leibnitz  Pule  on  the  pressure  gradient  term,  the  (BHS)  becomes, 


i  JL 

p  3x 


Pdz  + 


1  -h 


—  P 
P  n 


In 

3x 


l  _  3n  l 
—  P  — -  -f  —  x 
p  -h  3x  p  zx 


Assuming  that  the  pressure  at  the  free  surface  P^  is  zero  and 


realizing  that  x  and  x  are  the  time  averaged  surface  and 
zx  zx  . 

n  -h 

bottom  shear  stresses,  the  (RHS)  is  rewritten  as 


1  _3 
p  3x 


m 


Pdz  +  r-  P  ^ 


-h 


p  -h  3x 


+ 


1  - 

p  Tbx 


The  mean  pressure  at  the  bottom  P_h  can  be  defined  as  the 
sum  of  the  dynamic,  or  wave  induced  pressure  at  the  bottom,  and 
the  hydrostatic  pressure  at  the  bottom  such  that 

"  Pdyn  „  +  P5,h  *  *> 

-n 

_ 

Therefore,  P  —  can  be  written  as, 


—  in 

-h  3x 


+  pg(h  +  n) 


3h 

3x  ’ 


or  alternatively  as 
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Ip—  in 

p  -h  3x 


-  PJ 
p  dyn. 


3h  1  3  , 

9^r+  2  {g(h  + 


n) Z)  -  g(h  +  n) 


in 

3x  ‘ 


Using  the  result  the  time  averaged,  depth  integrated  x  momentum 
equation  is  given  by 


■srr  (u(h  +  n)  +  udz}  +  *r“  (u2(h  +  n)  +  2u  (  udz  +  [  u^dz 

3t  j_h  3x 


rn  -  r 

udz  + 

> -h  > — V 


2  p  ^  _  2  3  _ _  _ 

—  j  Pdz  -  y  g(h  +  n)  }  +  —  tuv(h  +  n)  +  U  |  vdz  +  v  I  udz 


3y 


-  r  »  -  rn 

U  vdz  +  v 
J-h  J-h 


rn 

+  uvdz} 
•'-h 


IP  il 

P  dyn  ,  3x 
— n 


-  g(h  +  n) 


in  +  I  — 

3x  p  sx 


—  T. 

p  bx 


u'  2dz 


3y  J  ,  u'v'dz 
— n 


(9) 


The  quantities  called  the  radiation  stresses,  or  the  excess 
momentum  fluxes  due  to  the  presence  of  waves  are  defined  as 
foil ows , 

S  =  [  (P+pu2) dz  -  ipg(h  +  r\)  2 - - -  {f  pudz}2 

J~h  p(h  +  n)  J-h 
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puvdz  - 


M" 

+  n) 


pvdz  pudz 


P (h  +  n)  J ~h 


(P+pv2)dz  -  ^  pg(h  +  n) 2  - 


p  (h  +  n)  J-h 


pvdz}' 


The  following  assumptions  will  also  be  used: 


(1)  The  product  of  the  wave  induced  pressure  at  the 
bottom  and  the  bottom  slope  will  be  assumed 
negligible, 

(2)  The  gradient  of  the  pressure  due  to  turbulent 

fluctuations ,  _ 

*  m 


a  r 

•*L“ 


'dz  ,  is  neglected,  and 


(3)  The  lateral  friction  caused  by  momentum  fluxes 
due  to  turbulent  fluctuations  is  assumed  to  be 
independent  of  depth. 


This  lateral  friction  will  be  called  x^  and  is  defined  as  x^  =  -pu'v1 
Finally  letting  the  total  depth  D  be  defined  as  D  =  (h  +  n) ,  the  x 
momentum  equation  in  its  final  form  can  be  written  as. 


at™  +  +  17  (UVD)  =  ~gD  £  •  ?-57 


,  3S  ,3S  ,  , 

1  — *Z  .  I  +  I  —  _  I  — 
p  3y  p  3x  p  sx  p  bx 


Following  the  same  procedure  the  final  form  of  the  y 


momentum  equation  can  be  found  to  be. 
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ay 


(VD) 


+  ta  (UVD) 


,  9  ,„2  .  9r| 

+  37  lv  Dl  ■  i; 


p  3x 


,9s  ,  as 

_  I .  -xx  .  1  __m  +  i  —  .  I  — 

P  3y  P  9y  p  sy  p  by 


(11) 


WAVE  TRANSFORMATIONS  INCLUDING  WAVE-CURRENT  INTERACTION 

The  equations  which  govern  both  wave  refraction  and  shoaling 
as  a  result  of  wave-current  interaction  used  in  the  model  are  those 
of  Noda,  et  al.  The  original  derivation  can  be  found  in  the  report 
presented  by  Noda,  et  al.  The  advantage  of  Noda's  method  is  that 
it  predicts  the  wave  angles  and  wave  heights  at  certain  points 
rather  than  along  a  wave  ray-  This  procedure  lends  itself  well  to 
use  in  the  finite  difference  model  because  calculations  are  performed 
at  points  which  lie  in  the  center  of  rectangular  grid  elements 
which  are  part  of  a  larger  grid  mesh. 

Starting  with  a  progressive  linear  gravity  wave,  the  free 
surface  can  be  written  as, 

n(x,y,t)  =  a(x,y,t) cos{<|> (x,y,t) } 

where  a  is  the  wave  amplitude  and  $  is  some  phase  function.  A 
wave  number  vector  can  be  defined  as 


k  =  V<|> 


(12) 


and  a  wave  scalar  frequency  can  be  defined  as 


a 


3 

at 


(13) 


Using  the  mathematical  property  that  the  curl  of  a  gradient  is 
identically  zero,  it  is  shown  that 


which  implies  that 


7  x  V<j>  =  0 


V  x  k  =  0 


This  equality  states  that  the  wave  number  vector  is  irrotational ;  i.e., 
the  wave  in  question  cannot  travel  in  circles.  Assuming  <j>  (x,y,t) 
is  continuous 

st  <7*>  -  7  ft 

or  substituting  Equations  (12)  and  (13)  into  the  above  expression, 
it  is  found  that 

- 

H  +  VC  =  0  (14) 

which  is  the  classical  conservation  of  waves  equation. 

For  the  case  of  a  wave  propagating  on  a  current  with  velocity 
\S  **  u^  +  vj,  it  earn  be  shown  that  the  scalar  frequency  with  respect 
to  a  stationary  reference  frame  is 


o  =  a  +  k  *  u 

The  wave  frequency  with  respect  to  a  moving  reference  frame  is 
given  by  the  dispersion  relation, 

a ^  =  gk  tanh  kh 

If  it  is  also  assumed  that  the  wave  number  field  changes  slowly 
with  time  then  from  Equation  (14) 

V(a  +  k  •  u)  =0 
or 

-f 

a  +  k  •  u  =  constant 

This  constant  can  be  evaluated  for  the  case  where  u  =  0  in  which 
2tt 

case  o  =  —  where  T  is  the  wave  period . 

Equation  (16)  then  becomes 


o  +  k  •  u  =  —  . 


Using  the  coordinate  system  shown  in  Figure  1  and  expanding 

Equations  (12)  and  (17)  into  Cartesian  coordinates  and  using  the 

dispersion  relation,  the  equations  which  govern  wave  refraction 

through  wave-current  interaction  are  given  by 

„  ,39  1  3k,  .  a  ,30  1  3k,  _ 

003  9  ^3x”  k  3y}  +  sln  9  {37+  k  3^}=  0 


{gk  tanh(kh)  }^2  +  uk  cos  0  +  vk  sin  0  *  . 


where  fl,k,h,u  and  v  are  all  functions  that  may  vary  in  both  the  x 
and  y  directions. 
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Figure  1  Planform  Definition  Sketch,  for  the  Wave  Transformation  Equations 

The  shoaling  of  waves  is  also  affected  by  the  interaction  of 
waves  and  currents.  The  effect  on  the  waves  is  determined  by  solving 
the  energy  equation.  The  form  of  the  energy  equation  used  in  the 
model  is  the  time  averaged  (over  one  wave  period) ,  depth  integrated 
form  for  the  case  of  a  wave  propagating  on  a  current  given  by 
u  =  ui  +  vj.  Turbulent  effects  and  dissipation  are  neglected  in 
the  derivation  which  can  be  found  in  Phillips  (1966) .  His  result 
is  given  by, 

(E,U+  cg*>>  +  ^7  tE'”  V1  +  Sxxf? 

+  S  |S+  S  |i  +  s  |r 

xy  3y  yy  3y  yx  3y 


0 


(20) 


Dividing  by  E  and  expanding  in  Cartesian  coordinates 


Equation  (20)  can  be  written, 


{■:!§■+  (u  +  C  cos  6)  7  |^  +  (v  +  C  sin9  ) 
E  3t  g  E  ox  g 


1  3E 
E  3y 


3  $ 

+  -r—  (u  +  C  cos  9)  +  t  (v  +  C  sin  0) 
dx  g  dy  g 


1  ,  3u  3u  9v  dv, 

+  —  (s  ~7. —  +  S  -5 —  +  S  —  +  S  -r — }  =  0 

E  xx  Sx  xy  3y  yy  3y  xy  3y 


In  the  above  equation,  E  is  the  total  energy,  both  potential  plus 
kinetic  of  a  progressive  linear  water  wave  and  is  given  by 

E  =  |pg  H2 

where  H  is  the  wave  height. 

Using  this  result,  carrying  out  the  differentiation,  and 
letting  a  quantity  Q  be  defined  as 

_  _  1  ,  3u  _  3u  „  3v  _  3v-i 

Q  =  ~  iS  -5 —  +  S  -5 —  +  S  -5 —  +  S  -z — } 

E  xx  dx  xy  dy  xy  dx  yy  dx 


the  energy  equation  becomes 
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For  all  applications  of  the  model  the  wave  height  H  is  assumed 

dH 

constant  in  time  so  —  =  0.  From  linear  wave  theory  the  group 

dt 

velocity  is  given  by 

=  C  ,  2kh  , 

g  2  sinh(2kh) t 

where 

C  =  tanh(kh)  }^//2 

is  the  wave  speed  or  celerity,  k  is  the  wave  number,  and  h  is  the 
water  depth. 


RADIATION  STRESSES 

In  the  derivation  of  the  momentum  equations  the  radiation 

stress  terms  S  ,  S  and  S  „  were  defined.  Those  forms  can  be 
xx  xy  yy 

simplified,  and  it  has  been  shown,  Longuet-Higgins  and  Stewart 
(1962) ,  that  to  second  order  for  a  progressive,  linear,  small 
amplitude  wave  they  can  be  approximated  by, 

Sxx  =  E{(2n-l/2)cos26  +  (n-1/2) sin20 } 

S  =  En  cos  0  sin  0 
xy 

Syy  ■  E{<2n-1/2>sin20  +  (n-l/2) cos20 } 

where  E  is  the  wave  energy,  0  is  the  wave  angle,  and  n  is  the 
ratio  of  group  velocity  to  wave  celerity: 
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where 


E  =  ^  pg  H2 

=  -3.  =  —  (1  +  - ^2 — ) 

C  2  '  sinh  2kh' 


k  =  wave  number  (=  — — ) 

L 

L  «  wave  length 
h  =  water  depth  and 

H  =  wave  height. 


WIND  STRESS 

The  capability  to  handle  a  wind  stress  was  retained  from  the 
work  of  Birkemeier  and  Dalrymple  (1975) .  The  wind  stress  is  included 
as  the  surface  stress  in  the  horizontal  momentum  equations.  No 
wind  effects  were  included  in  applications  of  the  present  version 
of  the  model.  A  brief  summary  of  the  wind  stress  formulation  is  given 
below. 


The  form  for  the  wind  stress  terms  was  assumed  to  be 


Tsx  "  pK'WlWx 


Tgy  *  P*|w|Wy 


A 
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where  w  is  the  magnitude  of  the  wind  speed  and  w  and  w  are  the 

x  y 

x  and  y  horizontal  components  of  the  wind  speed.  The  density  of 
water  is  p  and  the  wind  stress  coefficient  k  is  based  on  the  work  of 
Vein  Dorn  (1953)  and  assumed  to  be  a  function  of  the  wind  speed  such 
that 

tc  =  k,  for  w  <  w 
1  —  c 

w 

c  2 

K  =  K,  +  K  (1 - )  for  w  >  w 

1  2  w  —  c 

wc  is  a  critical  wind  speed  taken  as  14  knots  and  the 
coefficients  and  <2  are  taken  to  be  1.1  x  10  and  2.5  x  10  , 

respectively.  A  comparison  of  this  coefficient,  k,  with  real 
data  is  given  in  the  work  of  Pearce  (1972) . 


BOTTOM  STRESS 

The  bottom  friction  used  in  the  model  is  of  the  quadratic  form 


where  " - "  again  denotes  the  time  average  over  one  wave  period. 

The  total  velocity  vector  ut  is  comprised  of  both  mean  currents 
and  wave  orbital  velocities.  The  quantities  p  and  f  are  the  water 
density  and  the  Darcy-Weisbach  friction  factor,  respectively. 
Referring  to  Figure  2,  the  mean  current,  u,  can  be  broken  into  its 


x  and  y  components  u  and  v 


Figure  2  Definition  Sketch  for  the  Bottom  Friction  Formulation 


The  wave  angle  is  given  by  theta  and  1  and  ^  are  unit  vectors 
in  the  x  and  y  directions. 

The  total  velocity  vector  ut  can  be  expressed  as 
u^  =  (u  +  u  cos  0)i  +  (v  +  u  sin  0);) 
whose  magnitude  is  given  by, 

|ut|  =  \/u2  +  v2  +  u2  +  2uucos  0  +  2vusin  0 

A 

The  wave  orbital  velocity  u  can  be  expressed  as 

u  =  u  cos  o t  , 

m 

where  u  ,  the  maximum  wave  orbital  velocity,  is 
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m  T  sinh  kh 


The  x  and  y  components  of  the  time  averaged  bottom  friction 


then  become, 


t.  =  (  (u  +  u  cos0  cos  Jt)‘|u  |d(at) 

bx  16it  J  m  t 


_  f  r2ir  ,  ^ 

t  =  ■£—  (v  +  u  sinQ  cos  crt)  *  |u.  Id  (at) 
by  16ir  J  m  t' 


where  can  now  be  written  as. 


1^1  /  ^  2  a  2  2  A  a 

uj  =  vu  +  v  +  u  cos  at  +  2uu  cos  at  cos6  +  2vu  cos  at  sin0 
t1  m  m  m 


Both  of  these  stress  components  and  are  of  the  form 


— r 

16*  )r 


f  (ot)d(ot) 


Tb  "  llf  Sn 


where  S  is  a  sum  of  terms  to  approximate  the  integral  in  Equation  (24) 
n 


by  Simpson's  Rule: 


Sn  "  {fo(fft)  +  4fi(at>  +  2f2(at)  +  4f3(ot)  +  ... 


+  2f  (at)  +  4f  .(at)  -i  f  (at)} 
n-2  n-l  n 
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where  n  is  some  positive  even  integer  jnd 

Mat)  =  —  . 

n 

To  find  the  x  component  of  the  shear  stress  the  sum  can  be  found 
by  taking  the  integrand  in  Equation  (22)  for  the  function  f(ot). 
Similarly,  for  the  y  direction  the  function  f(at)  is  given  by  the 
integrand  in  Equation  (23)  .  The  value  of  n  was  chosen  to  be 
16.  This  choice  was  based  on  a  comparison  between  the  increased 
accuracy  achieved  with  an  n  value  greater  than  16  and  the  increased 
computation  time  associated  with  a  higher  n  value. 


LATERAL  MIXING 

In  a  paper  by  Longuet-Higgins  (1970) ,  the  author  presented 
a  formulation  for  the  steady  state  velocity  distribution  as  a 
result  of  waves  breaking  on  a  plane  beach  at  some  angle  to  the 
beach  normal.  This  formulation  was  based  on  a  balance  between  the 
bottom  friction  and  the  excess  momentum  flux  in  the  longshore 
direction  due  to  the  presence  of  waves.  The  resulting  longshore 
velocity  distribution  increased  linearly  from  zero  at  the  beach 
to  its  maximum  at  the  breaker  line.  Outside  the  breaker  line  the 
velocity  was  everywhere  zero. 

However,  physical  observation  and  both  laboratory  and  field 
data  indicate  that  mean  longshore  flows  are  present  beyond  the 


breaker  zone.  In  a  companion  paper  Longuet-Higgins  presented  a 
formulation  which  included  lateral  mixing  as  the  means  for  the 
alteration  of  the  linear  velocity  distribution  into  profiles  found 
in  nature  as  shown  in  Figure  3.  These  profiles  have  the  discontinuity 


Breaker  Line 
DISTANCE  OFFSHORE 

Figure  3  Longuet-Higgins*  Analytical  Solution  for  Oblique 
Wave  Attack  on  a  Plane  Beach 

at  the  breaker  line  removed  and  the  peak  velocity  shifted  shoreward. 
The  velocities  do  tend  to  zero  some  distance  outside  of  the  breaker 
line. 

There  is  a  physical  model  to  explain  the  mechanism  of  mixing 
and  it  is  based  on  momentum  exchange  between  fluid  elements  as 
they  fluctuate.  Consider  a  velocity  distribution  as  shown  in  Figure  4 
Allow  the  fluid  element  at  x^  with  mean  velocity  v^  to  have  a 
turbulent  fluctuation  u'  in  the  direction  parallel  to  the  x  axis. 


Figure  4  Definition  Sketch  for  the  Lateral  Mixing  Formulation 


If  this  fluctuation  u'  caused  the  fluid  element  to  move  to  its  new 
position  ,  then  the  element  would  be  accelerated  by  the  faster 
moving  fluid  with  velocity  V2 ,  thereby  increasing  the  momentum  of 
the  fluid  element.  The  flux  of  momentum  in  translating  fluid  from 
x^  to  x2  is  pu' .  Multiplying  by  -v'  which  is  the  difference  in 
velocities  between  the  two  fluid  layers  gives  the  momentum  change 
per  unit  time  in  the  direction  of  the  mean  flow  or  conversely,  the 
shear  stress  exerted  by  the  fluid  layer  at  x^  on  the  fluid  layer 
at  x2  given  by  t  *  -pu'v'.  The  negative  sign  is  a  consequence  of  a 
positive  turbulent  velocity  u'  causing  a  negative  shear  stress 
because  the  layer  at  x2  is  impeded  by  the  fluctuation  of  the  fluid 


element  from  layer  1. 
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Using  a  Taylor  series  approximation  to  first  order  between 
the  points  x^  and  Xj  to  find  an  expression  for  v'  gives. 


v' 


v2-vl 


(V1  +  *x  S*  -  V1  " 


SL 

x 


3v 

9x  ’ 


Therefore,  the  lateral  shear  stress  between  the  fluid  layers  can  be 
written  as, 


r  =  -pu ' 


l 


3v 
x  9x 


Since  for  an  arbitrary  coordinate  system  the  velocity 
distribution  could  vary  in  two  directions  the  shear  could  also 
include  a  term 


T  =  -pv' 


y  3y 


For  this  reason  the  lateral  shear  stress  x^  will  be  assumed  to  be. 


Two  coefficients  of  lateral  mixing,  one  for  each  direction  x  and  y, 
are  defined  such  that, 


e  5  u'Jt  and  e  =  v'  1 
x  x  y  y 


The  lateral  shear  stress  is  finally  written  as 


3u  dv 

T£=‘P  "y  97+ex^ 
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Longue t-Higgins  argued  that  the  mixing  coefficient  e^ 
must  tend  to  zero  as  the  shoreline  was  approached  since  the  turbulent 
eddies  responsible  for  mixing  cannot  be  of  a  scale  greater  than  the 
distance  to  shore.  He  assumed  that  e^  is  proportional  to  the  distance 
offshore,  x,  multiplied  by  some  velocity  scale  which  he  chose  to 
be  *^gh,  the  speed  of  a  wave  in  shallow  water  where  h  is  the  local 
water  depth.  Therefore,  e^  can  be  written  as 

E^  =  Nx  /gh 

where  N  is  a  dimensionless  constant  whose  limits  Longuet-Higgins  gave 
as 

0  <_  N  <_  0.016  . 

In  this  model  the  coefficient,  e^,  was  allowed  to  vary  linearly 
with  x  to  some  value  around  the  breaker  line.  From  this  point 
seaward  the  coefficient  remained  at  this  constant  value.  The  reason 
for  this  approximation  is  that  physically  there  must  be  some  limit 
on  the  scale  of  these  eddies.  This  limit  is  at  present  not  known. 

The  coefficient,  e^,  was  chosen  to  be  constant.  It  is  important  to 
note  that  the  purpose  of  the  model  is  to  present  a  stable  numerical 
scheme  which  includes  the  effects  of  mixing.  It  is  not  an  attempt 
to  verify  the  choice  of  the  mixing  coefficients  used. 


CHAPTER  III 
FINITE  DIFFERENCE 
FORMULATION 


At  this  point  we  have  established  a  set  of  three  governing 
equations,  both  mass  and  momentum,  to  be  solved  along  with  wave 
transformation  equations  and  forcing  terms  from  the  momentum  equations 
This  set  of  equations  cannot  be  solved  analytically  so  another  method 
of  solution  must  be  found.  In  this  chapter  a  means  for  solving  the 
mass,  momentum,  and  wave  trams formation  equations  numerically  is 
presented.  In  the  following  chapter  the  boundary  and  initial 
conditions  will  be  formulated  along  with  the  actual  solution  technique 

To  solve  the  problem  numerically  a  numerical  scheme  must  be 
defined  which  leads  to  a  systematic  method  of  solution.  This 
requires  the  following:  (1)  a  breakdown  of  the  area  under  study 
into  a  well  defined  grid  system  with  a  systematic  way  of  defining 
variables  of  interest,  and  (2)  the  conversion  of  the  governing 
equations  into  their  finite  differenced  forms. 

A  rectangular  grid  mesh  was  established  over  the  area  of 
interest  as  shown  in  Figure  5  where  x  and  y  denote  the  offshore  and 
longshore  directions,  respectively.  The  size  of  the  grid  blocks  is 
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j 


Figure  5  Grid  Schema  Definitions 

given  by  Ax  in  the  x  direction  and  Ay  in  the  y  direction.  At  each 

grid  block,  designated  by  the  subscripts  i  and  j  as  shown  in  Figure  6, 

certain  variables  must  be  defined.  The  velocities  u.  .  and  v,  . 

1/3  i/j 

are  defined  as  being  positive  if  they  enter  the  i,j  grid  block  in  the 


positive  x  and  y  directions.  Every  other  variable  will  be  defined 
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at  the  grid  center.  The  choice  of  this  definition  scheme  lends 
itself  well  to  the  boundary  conditions  used  in  the  model. 

The  methods  used  to  transform  the  continuity  and  momentum 
equations  into  their  differenced  forms  will  be  dealt  with  in  detail. 

The  derivation  for  the  other  equations  will  be  mentioned  only 
briefly.  Again,  the  more  detailed  derivation  can  be  found  in 
Birkemeier  and  Dalryraple.  Following  the  methods  of  Blumberg  (1977) 
and  Lilly  (1965) ,  certain  operators  are  used  throughout  the  derivation  of 
the  integrated  mass  and  momentum  equations.  The  first  two  operators  are 
essentially  central  finite  differences  and  are  given  by, 

6x(F(x,y,t)}  =  ^  {F(x+  ^-,y,t)  -  F(x-  ^j,y,t)} 

6x(F(x,y,t)}  =  {F(x+  Ax,y,t)  -  F(x-Ax,y,t)]  . 

In  the  first  case  the  differencing  takes  place  over  one  spacial 
grid  (Ax)  and  in  the  second  case,  over  two  full  grids  (2Ax) .  Also 
defined  are  the  following  two  operators  which  are  merely  averages 
in  space,  first  in  one  direction,  then  in  two. 

These  are  given  by, 

F(x,y,t)x  =  j  {F(x+  — ^-,y,t)  +  F(x-  ^,y,t) 


■v 


F(x,y,t) 


;y 


=  F(x,y,t) 
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Note  that  in  all  four  operators  ',  F(x,y,t)  is  some  arbitrary  function 
which  varies  in  space  and  time.  Similar  operators  also  exist  for 
time  t  and  in  the  other  horizontal  dimension,  y. 


The  governing  continuity  and  momentum  equations,  Equations  (6) , 
(10)  and  (11),  can  be  written  in  the  following  differenced  forms. 


CONTINUITY : 


6  (n  )  +  6  (D  u)  +  6  (Irv)  =  0 
t  x  y 


X  MOMENTUM: 


fit(DXu  )  +  6x(ixu  u*) 


6  (r?v  u7) 


-g  D*  <5  (n)  +  -  t  X  -  ~  t,  X  -  -  6  (i  Xy) 
x  psx  Pb*  py  xy 


-  p  VSxx>  +  ?  4y  SV“»>  »  ?  VS*  5x(v)1 


y  MOMENTUM: 


u  %  - y 

6  (D^v  )  +  6  (D*u  v^)  +  <5  (D^v  v^)  = 
t  X  v 


■g  7  6  (n)  +  —  t  y  -  —  t  ^  -  —  6  (S  ) 
y  P  sy  p  by  p  y  yy' 


?  Sx(5^*y»  *  ?  S«  ^  8x  <7*  V*” 


- v  - 


i  -  *  — 


(26) 


(27) 


(28) 


It  is  important  to  note  that  when  using  the  above  differenced 

equations  that  the  x,y  coordinate  is  defined  wherever  the  quantity 

to  be  solved  for  is  defined.  For  example,  since  the  x  momentum 

equation  is  used  to  solve  for  the  u  velocity  component,  the  x,y 

coordinate  is  located  where  u.  .is  defined.  In  the  continuity 

equation,  the  free  surface  displacement  n  is  to  be  solved  for  so 

the  x»y  coordinate  is  defined  to  be  at  the  center  of  the  grid 

where  n •  -is  defined.  This  is  important  in  converting  from  x,y 
1  /  D 

notation  into  i,j  notation. 


Figure  6  Differencing  Coordinate  (i,j)  Sketch 
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p  Tbx 


2p  ^Tbx  i,j+Tbx  i-1 , 


-  <S  (i  xy)  = 
p  y  xy 


<S.„.  ,  j., “S _  .  . +S, 


4Ayp  '°xy  i,j+l  ^xy  i,j-l'^xy  i-l,j+l-Sxy  i-l,j-l 


«  (S  ) 

X  XX 


.  ,-S 

XX  1,3 


XX 


i-1 


— x  „  r  ,  .  *  ,  1  r  n-1 . n-1  n-1,  n-1  n-1. 


n-1 .  n-1  n-1  . , 

e  (u.  .-u.  .  , )} 
y  i,j  i,3-l 


—x  .  f-xy,  1  ff  n— 1  n— 1  ■, , ,  n-1  n-1  n-1  n-1  . 

D  31  sAxAy  |^Di,j+Di-l,j}{  (ex  i,j+l  Ex  i,j  x  i-l,j  x  i-l,j+l) 

. .  n-1  n-1  .  ,  n-1  n-1  n-1  n-1  . .  ,  n-1  n-1  H 

*vi,j+l  vi-l,j+l  Ex  i,j  £x  i,j-l  £x  i-l,j-l  x  i-l,j  i,j  i-?»j  J 


y  MOMENTUM; 


-  ,ry  v  1  r  ,_n*l.  _n+l  .  n+1 


,  n-1  _n-l  .  n-li 
(D.  .+D.  .  .)v.  .} 

i,j  i,j-l 


-  is  EV)'°i.)iVi.)*Vj-iVi),»i,h]  • 
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{(v.  . . «+v.  . ) }-{ (D.  ,+D.  .  ,)V.  .+(D.  .  .+D.  .  ,)v.  .}’{(v.  .+v  .  ,)} 

i,j+l  i,j'  1,3  1,3-1  i*j-l  l»j"2  1,3-1  1*3  i,3“l  J 


-9  D*  «y(n)  »  -  jf7 


sy  i, 


.+t  .  .  . ) 

5  sy  i,3-l 


(31) 


—  7T* 

P  by 


2p  *Tbx  i, j+Tby  i,j-l* 


(S  )  = 

p  y  yy 


-j-  (s  .  .-s  .  .  .) 

pAy  yy  1,3  yy  1,3-1 


^  6  (F^Y)  = 

P  x  xy 


-  -  /c  _g  +g  —3  \ 

4pAx  xy  i+l,j-l  xy  xy  i+l,j  xy  i-l,j 


rv  r  ,  ,  ,  .  ,  1  ,_n-l. „n-l  .  r_  n-1,  n-1  n-1 

y  «x  Vy<“> 1  -  £5®i.J,l!i,H)ley  (ui*i,j-ui+i,j-i>  - 


n-1 .  n-1  n-1  . , 

e  (u.  ,-u.  .  ,)} 
y  1,3  i,j-l 


{e~*y6  (v)}«  — - — 5-  )  { (e  ?7i  .+e  **7,  H  ,+e  ?”*+e  ?"*_ 

xxx  8(Ax)2  ^  *  1+1'^  x  i+l»j-l  x  1»3  * 


. .  n-1  n-1. 


,  n-1 .  n-1  ,  n-1  ,  n-1  .  ,  n-1  n-1  .  71 

-(e  .  ,+e  .  ,  ,+e  .  ,  ,+e  .  ,  .  ,)(v.  .-v.  ,  .) } 

x  i,j  x  i, j-1  x  i-1, j  x  i-1, j— 1  1,3  1-1 »j  J 


The  superscripts  n,n+l,n-l  denote  the  time  level  of  a  particular 
quantity.  If  the  time  level  is  not  specified  it  is  assumed  to  be 
equal  to  n.  Also,  the  horizontal  mixing  terms  are  lagged  in  time 
for  stability  reasons,  as  mentioned  in  Holland  and  Lin  (1975). 
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Looking  back  at  the  set  of  equations.  Equations  (18) ,  (19)  and 

(21) ,  which  govern  the  refraction  and  shoaling  of  waves  through  wave- 

current  interaction,  if  Equation  (19)  is  differentiated  with  respect 
8k  8k 

to  x  to  get  and  with  respect  to  y  to  get  these  quantities  can 
be  substituted  into  Equation  (18) ,  which  can  then  be  written  in  the 
following  form: 


96  36 

A9^+B  97=C 


(32) 


where  A,  B  and  C  are  functions  of  the  quantities  sin  6,  cos  6,  k,  h,  u 

80 

and  v.  By  taking  a  forward  difference  in  x  to  approximate  and  a 

80 

backwards  difference  in  y  to  approximate  Equation  (32)  becomes: 


0.  . 

1/3 


D  +  E  0.  .  -  F  0.^,  . 

l,j-l  1+1,3 


(33) 


where  D,  E,  and  F  are  now  functions  of  the  quantities  sin  0,  ., 

ii] 

cos  0.  k,  h,  .,  u.  .,  v.  ..  To  evaluate  sin  9.  .  and 
i,3  i»j  i/3  i/3  i/3  i/j 

cos  j  Noda  used  a  first  order  Taylor  series  expansion  to  the  four 
neighboring  grids  i+l,j.  i-l,j,  i,j+l  and  i,j-l,  summed  the  results 
and  took  an  average  value. 

The  theta  field  0  is  solved  for  in  the  following  way. 

1/3 

Snell's  Law  is  used  to  approximate  the  angles  at  the  offshore  row. 
Working  shoreward  Equation  (33)  is  solved  for  in  a  row-by-row  relaxation 
until  the  angles  converge  to  their  correct  values  with  wave- current 
interaction  included.  After  each  updated  value  of  theta,  a  new  wave 
number  must  be  solved  for. 
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Equation  (19)  can  be  mitten  as 


,1/2 


2ir 


E (k)  =  {gk  tanh(kh) }  +  ukcos0  +  vksinQ  - —  =  0 


(34) 


To  solve  for  the  wave  number,  k,  after  each  updated  angle  is  found, 
the  Newton-Raphson  Method,  or  "method  of  tangents",  is  used.  This 
method  states  that 


k  E(w 

new  old  E*  (k  ) 
ola 


Differentiating  Equation  (34) ,  k  is  iteratively  solved  for  until 

Ik  -  k  , .001  |k  |  . 

1  new  old1  —  1  new1 


The  wave  height  field  is  calculated  in  much  the  same  way  as 

u 

the  wave  angle  field.  Multiplying  Equation  (21)  by  —  and  letting 
3H 

—  =0,  the  energy  equation  can  now  be  written  in  the  form 


.  3H  3H  _  „ 

A  — —  +  B  t —  *  C  H 

3x  3y 


(35) 


where  A,  B  and  C  are  functions  of  u,  v,  cos  8,  sin  8,  C  ,  Ax,  Ay  and 

9 

the  radiation  stresses.  Taking  a  forward  difference  in  x  to 

dH  3h 

approximate  ^  and  a  backward  difference  in  y  to  approximate  and 

solving  for  H .  ,  it  can  be  shown  that 

1rJ 


H.  .  -  D  H.  .  .  -  E 
i,j  i,j-l  i+1, j 


[ 

[ 


where  D  and  E  are  functions  of  the  same  quantities  as  A,  B  and  C. 

Again  the  offshore  row  of  wave  heights  are  obtained  from  shoaling 

and  refraction  due  to  Snell's  Law  and  the  wave  height  field  is 

determined  by  relaxing  row  by  row  in  the  shoreward  direction. 

On  each  row  a  solution  for  the  wave  height  is  reached  when 

|h  -  H  , J  <  .01  H  •  After  each  updated  value  of  H.  a 

1  new  old1  —  new 

breaking  wave  height  is  also  calculated  from  the  breaking  criteria 
given  by  the  Miche  formula 

(7O.  =  .12  tanh(kh)  . 

L  D  D 

If  the  calculated  H.  .is  larger  than  the  allowable  breaking  height, 

J 

the  height  H.  .  was  set  equal  to  the  breaking  height. 

1 1 3 


CHAPTER  XV 
METHOD  OF  SOLUTION 

In  the  previous  chapter  the  three  governing  equations, 
conservation  of  both  mass  and  horizontal  momentum,  were  derived  in 
their  finite  difference  forms.  The  numerical  procedure,  used  by 
Birkemeier  and  Dalrymple  and  Noda,  et  al. ,  to  calculate  the  wave 
characteristics  used  in  the  "forcing  terms"  in  the  momentum 
equations  was  also  presented.  In  this  chapter  the  initial  and 
boundary  conditions  will  be  stated  completing  the  problem  formula¬ 
tion,  whereupon  a  method  of  solution  will  be  given. 

In  every  application  of  the  present  model  the  initial 
conditions  were  assumed  to  be  the  state  of  rest.  The  velocity  field, 
both  u  and  v  components  as  well  as  the  mean  free  surface  displacement, 
n*  were  initially  set  equal  to  zero.  An  initial  depth  field  is 
inputted  into  the  model  and  the  wave  characteristics,  both  angles 
and  heights,  were  calculated  initially  using  Snell's  Law.  The  wave 
height  was  built  up  from  zero  to  its  full  deep  water  value  over  a 
certain  number  of  iterations  using  the  hyperbolic  tangent  function  in 
order  to  prevent  "shock  loading"  the  model.  The  form  for  this 
build-up  is. 
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H  =  H  tanh  (— ) 
o  T 

where 

Hq  =  deep  water  wave  height 
t  =  time  the  model  has  progressed 
T  =  time  allowed  for  wave  build-up. 

Certain  boundary  conditions  were  also  invoked.  Referring  to 
the  grid  system  shown  in  Figure  5,  no  offshore  flow  was  allowed 
into  row  m.  This  is  represented  by 


u 

m,] 


=  0 


which  essentially  simulates  a  wall  at  the  offshore  end  of  the  grid 
mesh.  The  reason  behind  this  choice  of  boundary  condition  is  this: 
at,  or  near,  steady  state,  if  row  m  is  far  enough  offshore  so  that 
the  effects  of  rip  currents  and  longshore  currents  are  negligible, 
the  assumption  of  zero  onshore-offshore  velocity  is  valid.  A  more 
correct  boundary  condition  would  be  to  let  the  mean  free  surface  be 
zero  in  deep  water  which  would  entail  the  addition  of  grid  blocks 
offshore,  thereby  increasing  computational  time  and  expense.  At  the 
inshore  end  of  the  grid  mesh  no  velocities  were  allowed  to  enter  the 
first  dry  grid  row.  Again,  this  implies  the  existance  of  a  wall  at 
the  beach  boundary. 


In  the  longshore  direction  periodic  boundary  conditions  were 


imposed.  Again  referring  to  Figure  5,  for  any  quantity,  Q,  periodicity 
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requires  that. 


Q(i,l)  =  Q (i ,N) 


Q(i,2)  =  Q(i,N+l) 

Q(i,3)  =  Q(i,N+2) 

and  so  forth.  Periodic  boundary  conditions  were  used  because  circu¬ 
lation  patterns  along  coastlines  often  tend  to  have  a  periodicity 
associated  with  them.  Also  if  it  is  desired  to  model  an  arbitrary 
stretch  of  beach  that  has  no  periodicity,  we  can  choose  boundaries 
far  enough  away  from  this  area  of  interest,  such  that  periodic 
boundary  conditions  in  the  longshore  direction  become  applicable. 


The  differenced  mass  and  momentum  equations  from  the  preceding 

chapter  were  derived  using  a  central  difference  in  time  for  the 

time  dependent  terms.  These  three  equations  can  also  be  written  in 
the  following  abbreviated  form. 


n+1  n-1  ,  _ 

n.  .  =  n.  .  +  2At  F. 
1,3  1.3  1 


(36) 


n+1 
u.  . 
i.3 


A  un_1  +  2At  F- 

1,3  2 


(37) 


n+1 


i-l 


v.  ;  =  b  v1?'  . 
1.3  1.3 


+  2At  F- 


(38) 


where  A  and  B  are  functions  of  the  depth  alone  and  F^,  F^,  and  F^ 
are  functions  of  all  the  variables  in  the  problem.  Also  remember  that 
Flf  F^.  and  F3  contain  quantities  evaluated  at  time  level,  n,  and  n-1. 
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These  three  equations  could  have  also  been  derived  using  a  forward 
difference  in  time  in  which  case, 


n 


n+l 
i,  j 


n+l 

u.  . 
i,3 

n+l 

v.  . 
i/D 


=  nn  .  +  it  f. 


=  C  u.  .  +  At  F 
1  ,  3 


=  D  v.  ,  +  At  F 
1,3 


4 


5  * 


(39) 


(40) 


(41) 


This  set  of  equations  will  be  used  to  start  the  computation  scheme 
as  will  be  seen  in  the  following  discussion. 

The  problem  has  now  been  formulated  with  a  set  of  three 
differenced,  governing  equations  (Equations  (36) ,  (37)  ,  and  (38) ) , 
initial  conditions  and  boundary  conditions.  The  next  step  is  to 
develop  a  computational  scheme  to  solve  the  problem.  Given  the 
initial  conditions,  the  wave  characteristics  and  the  values  of 
n,  u,  and  v  at  time  zero  (n*0) ,  the  functions  C,  D,  F^,  F^,  and  F^, 
in  the  set  of  Equations  (39)  through  (41)  ,  become  known.  Using  a 
time  step  of  l/2At  instead  of  At,  to  increase  the  accuracy  of  the 
first  calculations,  the  new  values  of  n,  u,  and  v  at  time  n+l  can 


be  calculated.  Knowing  these  new  values  the  wave-current  interaction 
equations  can  be  solved  to  find  wave  heights,  angles,  and  wave 
numbers  at  time  At/2. 


The  next  calculation  is  done  using  a  computational  scheme 
called  the  "leapfrog"  technique,  which  uses  the  set  of  Equations 
(36-38) .  To  employ  the  leapfrog  scheme  variables  must  be  defined 
for  the  two  previous  time  steps,  which  has  been  established  using 
the  initial  conditions  and  the  results  from  the  first  forward 
difference  calculation.  Referring  to  Figure  7,  the  variables  are 
known  at  time  levels  n-1  and  n.  The  functions  ,  F2,  F3,  A  and  B 
are  thus  known  and  u,  v,  and  n  at  time  level  n+1  can  be  calculated 
using  the  same  time  step  of  l/2At.  The  wave-current  equations  are 
again  solved  at  time  n+1  making  all  of  the  variables  defined  at  0 
and  At. 

time  level  n-1  n  n+1 

time  0  JjAt  At 


Leapfrog 


Figure  7  Initial  Forward  Difference  and  the  First  Leapfrog 
Solution  Steps 
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CHAPTER  V 
STABILITY  ANALYSIS 

The  exact  stability  criteria  for  the  full  sets  of  equations 
used  in  the  model  cannot  be  determined  analytically.  In  the  applica¬ 
tions  of  this  model/  stability  can  be  expressed  in  the  following  manner. 
The  speed  of  propagation  of  some  disturbance  in  the  model  must  be 
less  than  or  equal  to  the  speed  it  takes  the  disturbance  to  cross 
a  computational  grid  block  in  a  computational  time  step.  If  this 
criterion  is  not  met,  the  model  will  not  be  able  to  "see"  the 
disturbance . 

The  disturbance  speed  is  in  general  the  speed  of  a  shallow 
water  gravity  wave  plus  some  time  independent  mean  current. 

Therefore,  in  general,  the  stability  criteria  can  be  expressed  as, 

At<  V(Ax)2_+  (&y).2 _ 

|  VI  |  +  /gh 

The  maximum  magnitude  of  the  wave  speed  exceeds  the  minimum  current 
speed  in  general  so  the  stability  criterion  used  in  applications  of 


this  model  is, 
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(Ax)2  +  (Ay)2 


which  is  a  type  of  two-dimensional  Courant  number.  In  all  cases, 
the  time  step  actually  used  is  much  less  than  that  given  by  the 
above  criteria.  The  determination  of  a  better  stability  criteria 
and  the  probable  limits  on  the  time  step  were  not  investigated  in 
much  detail. 


As  mentioned  in  the  preceding  chapter  an  instability 
resulted  from  using  the  leapfrog  technique  to  integrate  the 
equations  in  time.  As  the  model  approached  a  steady  state,  the 
solution  diverged  into  two  disjoint  solutions;  one  associated  with 
the  even  time  steps  and  the  other  the  odd  steps.  These  solutions 
oscillated  with  growing  amplitudes  about  the  steady  state  solution. 
In  Roache  (1972) ,  the  author  referred  to  this  as  time  splitting. 


To  alleviate  the  problem,  a  leapfrog-backward  correction 
scheme,  Kurihara  (1965) ,  was  initiated  every  tenth  time  step.  The 
scheme  is  shown  below  as. 


h*  =  hn_1  +  2At  Gn 


hn+1  ■»  hn  +  At  G* 


where  h  may  be  u,  v,  or  rj.  Equation  (42)  is  simply  the  leapfrog 


calculation  done  by  the  Equations  (36-38)  where  *  denotes  the  new  or 


"interim"  time  level.  Using  the  new  values  u ,  v,  n  at  time  *,  the 
functions  G*  like  the  functions  F^,  T ^ ,  and  F3  from  Equations  (36-38) 
are  calculated  and  used  in  Equation  (43) ,  which  is  merely  a  backwards 
difference  in  time  to  the  level  n. 

This  scheme  was  chosen  because  it  damps  the  computational 
mode  of  the  solution  while  leaving  the  physical  mode  relatively 
unaffected.  For  a  more  in-depth  discussion  the  reader  is  referred 
to  the  work  by  Kurihara.  With  this  correction  scheme  included, 
which  essentially  "ties"  the  solutions  together,  every  tenth 
iteration,  the  model  proceeded  to  reach  a  steady  state  without  any 
further  time-splitting  instability. 


CHAPTER  VI 
RESULTS 


PLANE  BEACH  APPLICATIONS 

The  model  was  first  applied  to  the  case  of  a  single  progres¬ 
sive  wave  train  approaching  a  planar  beach  at  some  angle  to  the 
beach  normal.  Periodic  boundary  conditions  were  imposed  in  the 
longshore  direction  making  the  beach  infinitely  long.  The  purpose 
of  this  application  was  to  compare  the  model  to  the  earlier  efforts 
of  Birkemeier  and  Dalrymple  and  to  the  analytical  work  of  Longuet- 
Higgins  (1970)  on  longshore  currents  generated  by  obliquely  incident 
waves. 


For  the  series  of  plane  beach  runs,  the  same  input  data  was 
used.  The  deep  water  wave  characteristics  were:  (1)  wave  period 
of  8.0  seconds,  (2)  wave  angle  of  30.0  degrees  from  the  beach 
normal,  and  (3)  a  wave  height  of  2.0  meters  built  up  over  200 
(At  =0.5  second)  iterations.  The  time  step  of  0.5  seconds  is 
well  below  the  "allowable"  value  of  about  1.8  seconds  computed  from 
the  stability  criteria  presented  in  the  previous  chapter.  The 
region  of  interest  was  broken  into  a  6  x  30  grid  mesh  with  a  grid 
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spacing  of  10.0  meters  in  the  x  direction  and  15.0  meters  in  the 
y  direction.  The  beach  slope  was  chosen  to  be  0.025  which  resulted 
in  water  depths  of  0.0  to  7.00  meters.  In  all  runs  the  bottom 
friction  factor,  f,  was  chosen  to  be  0.08.  When  applicable,  the 
mixing  coefficients  N  and  were  chosen  to  be  0.01  and  0.5  meters/ 
sec2,  respectively. 


In  each  case  the  model  was  run  for  1200  iterations  which  is 
nearly  steady  state.  This  is  demonstrated  in  plots  of  mean  free 
surface  displacement,  or  velocity,  versus  time  at  selected  grid 
points  which,  for  the  case  using  the  present  model  without  horizontal 
mixing,  are  shown  in  Figures  9  through  12.  Note  the  dominant  period 
of  oscillation  of  about  112  seconds.  This  corresponds  well  to  the 
seiche  period  (of  113  seconds)  for  the  first  mode  of  oscillation  for 
a  triangular  basin  given  by  the  expression  (Wilson,  1966) , 


T  = 


3.28  L 


I'gh 


max 


where,  T  =  period  of  oscillation  in  the  basin 

L  =  length  of  the  basin 

h  =  maximum  depth  in  the  basin, 
max 


The  area  of  interest  for  the  plane  beach  applications  is  essentially 
an  infinitely  long  triangular  basin  as  a  result  of  the  flow  conditions 


imposed  at  the  inshore  and  offshore  boundaries. 
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Figure  10  Plot  of  Offshore  Meern  Free  Surface  Displacement  Versus  Time 

for  the  Present  Non-Linear  Model  Application  to  a  Plane  Beach 


f'8i  S3b  i(Tui3r§A*x  w 


Figure  11  Plot  of  Inshore  x-Velocity  Component  Versus  Time  for  the 
Present  Non-Linear  Model  Application  to  a  Plane  Beach 
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In  the  first  run  the  data  was  inputted  into  the  linear  model 
of  Birkemeier  and  Dalrymple.  The  resulting  longshore  current  profile 
is  shown  in  Figure  13.  Note  the  monotonic  increase  in  the  magnitude 
of  the  velocity  from  zero  at  the  dry  beach  to  a  maximum  of  about 
1.15  meters  per  second,  110  meters  offshore,  which  is  approximately 
the  location  of  the  breaker  line.  From  the  maximum,  the  velocity 
decreases  rapidly  to  zero  and  remains  essentially  zero  outside  the 
breaker  line,  indicating  no  offshore  mixing  effects. 

Note  the  similarity  between  the  form  of  this  current  profile 
and  the  form  given  by  Longuet-Higgins  (1970)  analytical  solution 
for  the  case  of  a  wave  approaching  a  planar  beach  at  some  oblique 
angle,  neglecting  the  effects  of  mixing,  as  shown  in  Figure  3.  The 
major  difference  between  the  two  profiles  is  that  the  linear  model 
result  shows  less  of  a  discontinuity  at  the  breaker  line.  This  is 
caused  by  the  use  of  a  discrete,  onshore-offshore  grid  size  in  the 
numerical  model  which  obscures  the  breaker  line.  As  this  grid 
dimension  is  reduced  the  location  of  the  breaker  line  becomes 
better  defined  thus  reducing  the  effect  of  "breaker  line  smoothing" 
on  the  velocity  distribution.  Also  associated  with  this  grid  size 
reduction  is  an  increase  in  the  magnitude  of  the  peak  velocity  to 
something  on  the  order  of  that  predicted  by  Longuet-Higgins,  which 
for  the  input  data,  is  about  1.5  meters  per  second. 


Figure  13  Longshore  Current  Profile  for  Plane  Beach  Case  from  Linear 
Model  of  Birkemeier  and  Dalrymple 
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A  lesser  point  to  note  is  that  in  all  the  plane  beach  applica¬ 
tions,  there  is  a  slight  discontinuity  in  the  velocity  profile  at  the 
inshore  region.  This  is  due  to  the  fact  that  initially  a  small  amount 

I 

of  water  was  placed  in  each  block  of  the  first  "dry"  grid  row.  As 
time  proceeded,  this  block  was  then  allowed  to  fill,  in  essence, 
simulating  the  effects  of  flooding.  Had  a  procedure  been  implemented 
that  would  have  allowed  for  the  flooding  of  more  than  one  grid  row, 
this  discontinuity  would  have  disappeared. 

Figure  14  shows  the  longshore  current  profile  resulting  from 
the  present  non-linear  version  of  the  model  excluding  the  effects  of 
horizontal  mixing.  The  form  is  similar  to  the  linear  result  except  for 
the  decrease  in  peak  velocity  and  the  extension  of  the  profile  a 
small  distance  outside  the  breaker  zone.  The  existence  of  velocities 
outside  the  breaker  zone  is  due  to  the  fact  that  the  advective 
acceleration  terms  in  the  differenced  y  momentum  equation  caused 
velocities  outside  the  breaker  zone  to  be  calculated  using  velocities 
inside  the  breaker  zone.  The  deviation  between  peak  velocities 
predicted  by  the  two  models  is  attributed  to  the  use  of  the  quadratic, 

"exact"  bottom  friction  formulation  used  in  the  non-linear  model. 

Liu  and  Dalrymple  (1978)  showed  that  the  use  of  this  bottom  friction 
term  caused  the  peak  velocities  of  Longuet-Higgins  linear  approximation 
to  be  decreased  by  about  20%  for  a  breaking  wave  angle  of  around 


10  degrees.  In  the  plane  beach  applications  of  the  linear  and  non¬ 
linear  models  the  breaker  angle  is  about  12  degrees  and  the  discrepancy 


Figure  14  Longshore  Current  l/roirle  for  Plane  Beach  Case  from  Present 
Model  Excluding  Mixing 


62 


in  peak  velocities  between  the  two  models  is  about  that  same  20%. 

The  final  run  using  the  plane  beach  data  was  made  with  the 
non-linear  model  including  lateral  mixing.  The  steady  state  long¬ 
shore  velocity  profile  is  shown  in  Figure  15.  This  profile  exhibits 
the  following  differences  from  that  predicted  by  the  model  neglecting 
mixing : 

(1)  in  the  inner  one  half  of  the  surf  zone  the 
velocities  are  increased  slightly, 

(2)  the  peak  velocity  is  shifted  to  a  new 
location  shoreward  of  the  breaker  line,  and 

(3)  the  velocity  distribution  extends  well 
beyond  the  breaker  line  eventually  decreasing 
towards  zero. 

The  form  of  this  profile  is  similar  to  the  analytic  result  of 
Longuet-Higgins  including  the  effects  of  mixing  shown  in  Figure  3. 

BARRED  PROFILE  APPLICATION 

Since,  in  nature,  beach  topography  is  often  comprised  of 
fragmented  longshore  bars,  the  present  version  of  the  model  was 
run  on  a  bottom  configuration  that  included  an  infinitely  long 
longshore  bar.  A  comparison  of  the  barred  profile  with  a  planar 
beach  (with  a  slope  of  0.025)  is  shown  in  Figure  16.  The  remaining 
input  into  the  model  was  identical  to  that  used  in  the  series  of 
plane  beach  runs.  The  model  was  run  both  with  and  without  the  effects 
of  mixing  included. 


The  results  for  the  case  without  mixing  are  shown  in 
Figure  17.  Notice  the  two  distinct  regions  where  a  longshore 
current  distribution  exists.  The  velocity  "spike"  offshore  is  due 
to  waves  breaking  on  the  bar.  As  the  wave  height  decreases,  as  a 
result  of  breaking,  an  onshore-offshore  gradient  of  y  momentum  flux 
is  created  which  drives  a  longshore  current.  In  the  trough,  however, 
the  wave  height  starts  to  reform  (no  more  breaking)  resulting  in 
the  absence  ^f  a  longshore  current  in  this  region.  In  reality,  a 
longshore  current  does  exist  in  the  trough,  Allender,  et  al.  (1978) , 
due  to  the  mechanisms  of  turbulent  dissipation  during  breaking 
within  a  bore,  lateral  mixing  which  has  been  included  in  the  model, 
and  a  set-up  of  water  within  the  trough,  Dalrymple  (1978) . 

Figure  18  shows  the  resulting  longshore  current  profile  for 
the  model  run  including  horizontal  mixing.  The  effects  of  mixing 
are  very  much  evident  as  the  amplitude  of  the  longshore  velocity 
"spike"  is  reduced,  the  discontinuities  in  the  velocity  profile 
disappear,  and  a  longshore  current  now  exists  in  the  bar  trough.  Had 
the  turbulent  energy  dissipation  mechamism  been  included  in  the 
model,  the  results  would  probably  have  approached  those  found  in 
nature . 

PERIODIC  BOTTOM  TOPOGRAPHY  APPLICATION 

The  model  was  next  applied  to  the  periodic  bottom  topography 
developed  by  Noda,  et  al.,  which  is  essentially  a  channel  at  some 


art 


angle  to  the  beach  normal.  The  formulation  for  this  bottom 


configuration  is  given  in  Appendix  A.  The  present  version  of  the 
model,  including  the  effects  of  mixing,  was  compared  to  the  linear 
model  of  Birkemeier  and  Dalrymple .  The  following  wave  characteristics 
were  used  in  both  instances: 


(1)  deep  water  wave  height  of  1.0  meters, 

(2)  wave  period  of  4.0  seconds,  and 

(3)  a  deep  water  wave  angle  of  30.0  degrees  to 
the  beach  normal . 


The  bottom  friction  factor,  f,  was  chosen  to  be  0.08,  and  the  mixing 

2 

coefficients,  N  and  e^,  were  chosen  to  be  0.005  and  0.5  m/sec  ,  respectively. 
In  both  runs  the  wave  height  was  built  up  to  its  deep  water  value 
over  100  seconds. 


Both  models  were  run  until  they  reached  approximately  a 
steady  state,  about  500  seconds.  The  wave-current  interaction  process 
was  halted  in  the  linear  model  after  150  seconds  because  the 
offshore  velocity  components  grew  too  large  for  the  refraction 
routines  to  handle.  In  the  non-linear  model,  however,  the  wave- 
current  interaction  process  was  included  for  the  duration  of  the 
run  time.  The  circulation  patterns  after  500  seconds  sure  shown  in 
Figures  19  and  20. 

Note  the  strength  of  the  rip  and  its  offshore  extent  in 
the  linear  model  compared  to  non-linear  model.  The  peak  velocity  in 


Figure  19 


Current  Vector  Plot  for  the  Model  of  Birkemeier  and 
Dalrymple  Run  on  the  Periodic  Bottom  Topography 
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Figure  20  Current  Vector  Plot  for  the  Present  Model  Including 
Mixing  Run  on  the  Periodic  Bottom  Topography.  The 
arbitrary  longshore  mixing  coefficient  probably  chosen 
as  too  large. 
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the  linear  model  run  is  about  3.0  meters  per  second  where,  as  in  the 
non-linear  model,  it  is  about  0.8  meters  per  second.  This  large 
discrepancy  is  due  to  the  inclusion  of  mixing  in  the  non-linear 
model.  The  mixing  tends  to  spread  the  rip  out  and  decrease  its 
offshore  velocity  components  thus  causing  the  rip  to  turn  more  in 
the  longshore  direction  as  shown  in  Figure  20.  The  effects  of  the 
convective  acceleration  terms  are  not  clearly  visible  because  it 
seems  as  though  the  form  of  the  rip  itself  is  governed  primarily  by 
horizontal  mixing. 

INTERSECTING  waves  application 

The  final  application  of  the  model  was  to  the  case  of  inter¬ 
secting  wave  trains  of  the  same  frequency  on  a  plane  beach  which 
Dalrymple  (1975)  showed  could  generate  rip  currents.  The  purpose 
of  this  application  was  to  investigate  the  effect  of  the  convective 
acceleration  terms  in  the  model.  The  following  derivation  closely 
follows  the  work  of  Dalrymple. 

Given  two  intersecting  wave  trains  A  and  B  with  amplitudes 
a  and  b  and  a  common  frequency,  a,  in  terms  of  the  coordinate  system 
shown  in  Figure  21,  the  free  surface  displacements  for  the  two 
wave  trains  can  be  written  as, 

*  a  cos(k  cos  ax  +  k  sin  ay  +  at) 

n2  * 


b  cos(k  cos  6x  +  k  sin  fJy  +  at) 


Figure  21  Definition  Sketch  for  the  Intersecting  Waves  Application 


The  total  free  surface  +  n2  can  then  be  written  as, 

nT  =  2a  cos  {^-(cosa  +  cos  g)x  +  -^(sina  +  sin  g)y  +  at} 

k  k 

cos  {ytcosa  -  COS  g)x  +  y(sina  -  sin  g)y} 

+  (b-a)cos  {k  cos  gx  +  k  sin  gy  +  at}  .  ( 

Using  the  linearized  dynamic  free  surface  boundary  condition  the 
velocity  potential  $T  can  be  shown  to  equal, 

$  =  -  sin{y(cosu  +  cos  g)x  +  y(sino  +  sin  g)+  at} 

T  a  cosn  Kn  i 

k  k 

cos  {^(coso  -  cos  g)x  +  j(sina  -  sin  g)y} 

+  cosh  k(h+z)  CQ8  gx  +  ^  8in  gy  +  0t} 
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Prom  the  velocity  potential  the  total  orbital  velocities  can  be 
found  from, 


3$ 

U  =  ”  3x 


_  3£ 

3y 


3^_ 

3z 


The  radiation  stresses,  which  are  essentially  the  forcing  terms,  are 
defined  as. 
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where 


P  =  pg(n 


-  3  "  -  2 

puw  dz  +  —  pvw  dz  -  pw 
z  ^  '  z 


Through  tedious  calculations  the  radiation  stresses  are  found  to  be, 

Sxx  "  4sinh2kh  [®2cos2ot  +  b2cos26  +  2ab  cosocosfcos{2(T)}J‘  {2kh+sinh2kh} 


-  (cos8-cosa) 2cos{2(5) }’ {2khcosh2kh-sinh2kh} 


-  8sinh2kh  (sina-sing) 2cos{2  >’ {2khcosh2kh-sinh2kh> 
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4sinh2kh~  (j*2+b2+2ab  cos{2(2)  }J- {sinh2kh-2kh) 


+  pgabcos2  (2)  +  7-  pg(b-a)' 


yy  4sinh2kh 


inh2kh  j^2s*-n2a+b2s*n28+  2absinasin6cos{2^2)}J‘ {2kh+sinh2kh} 


-  gg?^2kh  ^cos^"cosa^  2cos{2^2)  }*  {2khcosh2kh-sinh2kh} 


-  (sina-sin$) 2cos  (2(2)}*  {2khcosh2kh-sinh2kh} 

8sxnh2kh  w 


4sl^2kh  G2+b2+2abCOS  *2@  }3‘^sinh2kh_2kh} 


+  pgabcos  (2)  +  pg(b-a) 

inh2kh  j^2sinacosa+b2cos^sin^+abcos^2^2)  J  sin(.a+gTJ 


xy  4sinh2kh 


{ 2kh+sinh2kh} 


where  the  expression  "(jt)"  is  defined  as, 


© 


Ir  U 

2)  =  -jfcosa-cosPJx  +  j(sina-sinB) y 


The  tine  independent  mean  free  surface  displacement,  n, 
is  defined  by 


where  " - "  denotes  the  time  average  over  one  wave  period. 

Substituting  the  expressions  for  the  velocity  components  u,  v,  and 
w  from  the  velocity  potential  4>T,  n  can  be  written  as. 


n  = 


-k 


2sinh2kh 


|a^+b2+2abcos{2  (2) }  *  (cos  (a-(3)  cosh^kh-sinh2khJ  (45) 


where  "(2)"  is  the  same  quantity  defined  previously.  Notice  that 
the  mean  free  surface  displacement  is  modulated  in  the  x  and  y 
directions  by, 

cos{k(cosa-cos(3)  x  +  k(sina-sin$)y} 

Using  Snell's  Law  which  says, 

k  sina  =  k  sin  a  and  k  sin$  =  k  sin  B 
00  00 

and  using  the  fact  that  k  =  7—  ,  where  "o"  denotes  deep  water  values 

O  Li 

O 

for  the  wave  length,  L,  and  the  wave  angles,  a  and  B,  we  see  that 


there  is  a  periodicity  of  the  mean  displacement  in  the  longshore 
direction  with  a  periodic  spacing,  l,  given  by. 
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This  periodicity  in  water  level  and  wave  height  causes  water  to 
be  driven  from  regions  of  high  mean  water  level  displacement  to 
regions  of  lower  displacement  resulting  in  the  formation  of 
circulation  cells. 

In  order  to  attempt  to  model  this  phenomena,  certain 
simplifications  to  the  model  had  to  be  made.  Since  the  refraction 
and  shoaling  routines  borrowed  from  the  work  of  Noda,  et  al . ,  could 
not  treat  more  than  one  wave  train,  they  were  replaced  with 
routines  governed  by  Snell's  Law  neglecting  wave-current  interaction. 
Again  a  quadratic,  "exact,"  bottom  friction  was  used  including 
velocities  due  to  both  mean  currents  and,  this  time,  the  two  wave 
trains.  In  the  momentum  equations  the  advective  acceleration  terms 
were  retained,  horizontal  mixing  was  neglected,  and  the  radiation 
stresses  were  calculated  using  the  results  presented  earlier  in  this 
section.  A  procedure  for  determining  the  wave  heights  for  use  in 
the  radiation  stresses  and  the  bottom  frictional  stresses  is  given 
in  Appendix  B. 

Three  runs  were  made  using  different  combinations  of  wave 
heights  and  wave  angles.  The  remainder  of  the  input  data  for  all 
three  runs,  however,  was  the  same  and  is  given  as  follows.  The 
waves  were  run  on  a  plane  beach  with  a  slope  of  0.025.  The  planform 
area  of  interest  was  comprised  of  25  grids  in  the  x  direction  with 
an  Ax  grid  size  of  5.0  meters,  and  21  grids  in  the  longshore  direction 
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with  a  Ay  grid  size  of  4.0  meters.  The  time  step  was  chosen  to  be 
0.2  seconds  and  the  model  was  run  for  1500  iterations  for  all  three 
cases.  A  wave  period  of  7.159366  seconds  was  used  which  resulted 
in  rip  spacings  of  80.0  meters.  The  bottom  friction  factor  was 
set  equal  to  0.12  to  allow  the  system  to  reach  steady  state  after  the 
1500  iterations  and  to  decrease  the  magnitude  of  the  resultant  currents. 

The  first  case  run  used  waves  of  equal  heights  and  equal 
angles  on  either  side  of  the  beach  normal.  The  deep  water  wave 
heights  were  0.25  meters  and  the  deep  water  angles  were  +30.0 
degrees.  For  this  case,  referring  to  Equation  (44) ,  a=b  and  a=-B 
resulting  in  a  free  surface  displacement  given  by 

n  =  2a  cos(k  sin  a  y)cos(k  cos  ax  +  at) 

This  free  surface  describes  a  wave  train  moving  in  the 
-x  direction  with  a  modulated  wave  height  that  is  periodic  in  the 
longshore  direction  only.  The  periodicity  in  wave  height  is  the 
driving  mechanism  producing  the  rip  current  perpendicular  to 
the  beach  as  shown  in  Figure  22.  Note  the  constricted  width  of  the 
rip  current  in  relation  to  the  width  of  the  inflow  region.  This 
is  a  result  of  the  convective  acceleration  terms.  Also  note  the 
weak  rip  head  where  the  currents  diverge  from  the  rip  axis  and 
return  towards  shore. 


Figure  22  Current  Vector  Plot  for  a  Rip  Current  Perpendicular  to  the 
Shoreline 


In  the  second  case  the  waves  again  were  chosen  to  be  of 
equal  height  (0.25  meters)  but  one  wave  approached  normal  to  shore 
while  the  other  approached  at  an  angle  of  60.0  degrees  to  the  beach 
normal;  i.e.,  a=b  but  <*/-&.  The  resulting  free  surface  is  modulated 
in  both  the  x  and  y  directions  which  should  result  in  a  rip  current 
at  an  angle  to  the  beach  normal  given  by, 

«  .  -1  ,sin  a  -  sin  6. 

6  =  tan  ( - £■) 

cos  a  -  cos  0 

For  this  case  the  angle  should  be  about  30.0  degrees.  The  results 
are  shown  in  Figure  23. 

In  the  third  case  the  waves,  A  and  B,  had  different  heights, 
0.1  and  0.4  meters,  and  wave  angles  of  30.0  and  -30.0  degrees, 
respectively.  The  resulting  circulation  pattern  is  shown  in 
Figure  24.  Note  the  meandering  current  with  alternating  regions 
of  strong  and  weak  longshore  velocity  along  the  beach.  This 
circulation  would  lend  itself  well  to  the  formation  of  rythmic 
beach  features.  Looking  at  Equation  (44) ,  we  see  that  there  is  a 
non-zero  term, 

(b-a)cos(k  cos  6x  +  k  sin  6y  +  at} 

which  is  a  wave  train  at  an  angle  to  the  beach  normal  with  height 
2(b-a).  This  wave  is  present  in  addition  to  the  normal  wave  train 
with  the  modulated  height  from  the  first  case  causing  a  longshore 
current  which  is  superimposed  on  t) a  cellular  circulation. 


79 


Figure  23  Current  Vector  Plot  for  a  Rip  Current  at  an  Angle  to 
the  Shoreline.  Note  that  in  this  case  the  offshore 
wall  has  influenced  the  flow. 


0YM4.O  H 
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CHAPTER  VII 
CONCLUSIONS 

A  model  that  can  accurately  predict  currents  and  wave 
transformations  in  the  nearshore  zone  is  a  necessary  step  in  attempt¬ 
ing  to  predict  actual  changes  to  our  coastlines.  From  the  work 
done  in  this  pro ject  and  the  results  found,  it  appears  that  the 
inclusion  of  the  convective  acceleration  terms  and  lateral  mixing 
terms  in  the  horizontal  momentum  equations  have  important  effects 
on  models  used  to  predict  nearshore  circulation.  The  terms  become 
especially  significant  in  attempts  to  model  circulation  over 
irregular  bottom  topographies  which  include  bars  and  channels. 


There  are  still  many  aspects  of  the  model  that  could  be 
changed  to  make  it  even  more  complete.  Among  them  ares 


(1)  enable  the  model  to  handle  more  than  one  wave 
train  within  the  wave-current  interaction  process, 

(2)  include  a  better  wave-current  interaction  scheme, 
especially  in  the  surf  zone,  which  could  treat 
strong  offshore  flows  that  oppose  the  wave 
orbital  velocities, 

(3)  treat  the  continuity  and  momentum  equations  with 
an  implicit  scheme  to  avoid  stability  problems ,  or 
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(4)  obtain  a  better  knowledge  of  the  stability  of 
the  present  model  as  a  function  of  various 
parameters,  and 

(5)  create  the  ability  to  handle  a  variable  grid 
size  so  that  special  regions  of  interest  could 
be  investigated  in  more  detail. 


In  its  present  form  the  model  is  a  very  powerful  tool  yet  continual 
improvements  are  needed  to  maintain  its  usefulness. 


APPENDIX  A 

PERIODIC  BOTTOM  TOPOGRAPHY 

The  periodic  bottom  profile  used  in  the  model  was  developed 
by  Noda,  et  al.  (1974) .  The  depths  are  given  by 

h  =  rax  j\+A  exp{-3(-^)1//3}  sin10  ^  (y-x  tan  3)J 

where  m  =  beach  slope  =  .025 

x,y  are  the  coordinates  of  the  depth  location 
X  =  length  of  periodic  beach  =  80  meters 
A  =  amplitude  of  bottom  variation  »  20 
6  =  angle  of  rip  channel  to  beach  normal  =  30  degrees. 

The  grid  spacing  was  chosen  to  5.0  and  4.0  meters  in  the  x  and  y 
directions  respectively.  There  were  25  and  21  grids  in  the  x  and  y 
directions.  The  last  grid  row  and  the  "dry"  grid  rows  were  made  planar 
with  the  slope  being  .025. 


APPENDIX  B 


DETERMINATION  OF  WAVE  HEIGHTS  FOR  THE  CASE  OF  INTERSECTING  WAVES 

Given  the  two  wave  trains 

=  a  cos(k  cos  ax  +  k  sin  ay  +  at) 
n2  =  b  cos(k  cos  0x  +  k  sin  8y  +  at) 

The  total  free  surface  nT  =  +  n2  can  be  written  as, 

nT  =  a  cos(e+at)  +  bcos(6+ot) 

where  e  and  6  are  phase  functions.  Expanding  the  expression  for 
nT  and  rearranging, 

nT  =  cos  at  {a  cos  e+  b  cosfi)  -  sin  at  {a  sin  e+  b  sin  5} 


Defining 


and 


n_  becomes 
T 


D  cos  0  =a  cos  e  +  b  cos  6 


D  sin  6  =a  sin  e  +  b  sin  6 


nT  -  D  cos (ot  +  0) 


84 


85 


where 


D  =  v  (a  cos  e  +  b  cos  6)^  +  (a  sin  e  +  b  sin  6)^ 


Therefore,  the  total  wave  height  due  to  the  intersecting  waves  is 
2D.  Substituting  in  the  expressions  for  e  and  6,  the  total  wave 
height  can  be  written, 

H_  =  2.0 '/a2+b2+2ab  cos (k(cosa-cosB)x  +  k(sina-sin8)y3 


The  amplitudes  a  and  b  are  given  by  the  expressions 


a  =  —  H  k  k 
2  *  rA  S» 


b  =  —  H  ie  ic. 

2  B  rB  SB 


where  Kr  and  *s  are  the  refraction  and  shoaling  coefficients  derived 
using  Snell’s  Law. 


To  determine  if  breaking  occurs,  this  total  wave  height  is 
compared  to  a  limiting  breaker  height  of  <h,  here  k.  is  a  breaking 
coefficient  chosen  as  0.78.  If  the  total  wave  height  exceeded  the 
value  given  by  the  breaking  criterion,  that  value  of  <h  was  substi¬ 
tuted  into  Equation  (  46)  for  H^.  A  parameter  6  was  defined  such 
that  b  =  3a.  This  result  was  substituted  into  Equation  ( 46  )  for 
b.  The  equation  was  then  solved  for  a  which  gave  b.  The  wave 


heights  in  the  surf  zone  were  then  found  using  Equations  (47)  and  (48) 
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